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ABSTRACT 

Planetary systems form in gas-dust protoplanetary discs via the growth of solid 
bodies. In this paper, we show that the most intriguing stage of such growth — namely, 
the transformation of 1-10 m boulders into kilometre-sized planetesimals — can be 
explained by a mechanism of gravitational instability. The present work focused on 
the origin of self-gravitating clumps in which planetesimal formation could take place. 
Our computer simulations demonstrated that such clumps of gas and boulders formed 
due to the development of a two-phase instability. This instability revealed a so-called 
'mutual influence effect' in the protoplanetary disc, where the dynamics of the system 
were determined by the collisionless collective motion of a low-mass subdisc composed 
of primary solids. We found that a 0.1 Cs velocity dispersion in the boulder subdisc 
was sufficient to cause the formation of self-gravitating clumps of gas and boulders. 
In such regimes, the time needed for the formation of the collapsing objects was less 
than the boulders' dissipation time in the density waves of the medium. 
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1 INTRODUCTION 

Results from numerical modelling and observations of the 
gas-dust discs around young solar-type stars provide the ba- 
sis for the currently accepted scenario for the formation of 
the Solar System (Hartmann 2009). This model suggests 
that the disc forms simultaneously with a protostar, follow- 
ing the gravitational collapse of gas in the molecular cloud 
(Petit & Morbidelli 2005). In the first stage (which lasts up 
to 100,000 years), a protostar forms; this consists mainly of 
hydrogen and helium, and has a mass of about one tenth of 
the solar mass. 

The massive gas-dust disc is formed by the collision 
of the opposing gas streams. The gas streaming during 
the molecular cloud collapse is supersonic (Spitzer 1978). 
The infalling gas streams collide and produce the diverging 
shock waves (Landau & Livshitz 1987)that decelerate the 
gas streams velocity (see Fig. la). Between a pair of shock 
waves the gas density is higher than its extreme value for 
a single shock wave, as it was shown theoretically and ex- 
perimentally by Orishich, Ponomarenko & Snytnikov 1989. 
When the flow rate of collapsing gas decreases in the end 
of the star formation stage, the shock waves diverge from 
the disc plane (see Fig. lb). Now the gas can leave the disc. 
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Fast diminution of its flow provides conditions for the fast 
gas expansion which is followed up with its cooling. 

The molecular cloud dust moves together with the gas. 
During the gas compression behind the shock waves (Fig. la) 
the dust grains grow in size as a result of heavy molecular 
absorption and coagulation (Cuzzi & Weidenschilling 2006). 
Adsorption time is less than 10^ years in a medium with 
the hydrogen concentration higher than 10^ cm~^ (Spitzer 
1978). In the disc the concentration of gas is by several or- 
ders of magnitude higher than 10^ cm~^. So the molecules 
are adsorbed according to their sublimation temperature. 
Spitzer (1978) estimated that the coagulation length exceed 
the absorption length in the order of 1-2. These lengths are 
signiflcantly less than the massive disc thickness. That is 
why the dust grains can grow in size in accretion disc, set- 
tle on the disc plane and form their own subdisc (Cuzzi 
& Weidenschilling 2006). Due to the subdisc gas expansion 
in outside flow the relation of solid density to gas density 
increases in contrast to the molecular clouds, where dust 
constitutes 1-2 per cent of the mass. 

The protostar increases its mass up to the star mass due 
to accretion from the disc and the residues of the surround- 
ing molecular cloud. Alteration of the disc mass is deflned 
by the flows: (a) falling to the disc from the cloud, (b) ac- 
creting from the disc to the protostar, (c) leaving the disc, 
e.g. in jets. By the time discs become observable (1-3 million 




years) , their masses have decreased to 1 per cent of the host 
star mass. Inside the period between protostar and young 
star a stage of massive disc exists, when the mass of the disc 
and the mass of the central body are comparable. Over the 
course of the next 60-100 million years, the circumsolar disc 
evolves to a state similar to that currently present in our 
own Solar System. 

The formation of planets in circumstellar discs proceeds 
via the growth of solids; nanometer-sized dust in molecular 
clouds transforms into planets with radii of some thousands 
of kilometres. 

The growth of dust grains via coagulation depends on 
chemical composition. Usually ice, silicates and iron are 
mentioned as dust compounds. But due to the cosmic abun- 
dance of chemical elements (Lodders & Fegley 1998) the 
main components of disc solids dust should be organic com- 
pounds from H, C, O, N with high molecular weights (Herbst 
& van Dishoeck 2009). Such compounds can be easily syn- 
thesized on Mg, Si, O, Fe dust grains that work as catalysts 
in a medium with abundance of H and He (Khassin & Snyt- 
nikov 2005). That is why conglutination of multi-component 
organic compounds and non-organic bodies can be expected 
up to 1-10 meter in size. When the subdisc consists of such 
bodies with velocities greater than 10 km s~^ they collide 
less than once per orbital time. 

The growth of planetesimals (i.e., bodies with a size 
in excess of a kilometre) may occur due to collisional ac- 
cumulation, since their own gravitational field can attract 
smaller bodies and hold these to their surface (Safronov 
1969). The mechanism by which kilometre-sized planetesi- 
mals are formed from 'boulders' (or agglomerates, with sizes 
from 1 — 10 m) has been the subject of intense interest and 
discussion, and has formed the basis of many studies (Youdin 
& Shu 2002, Makalkin k Zighna 2004, Rice et al. 2006, Jo- 
hansen et al. 2007). At this stage, the concentration of solids 
provides motion, with rare collisions at orbital time with ve- 
locities of tens kilometers per second. The metre-sized boul- 
ders are not enlarged during high- velocity collisions. Colli- 
sions of such boulders don't result in sticking irrespective of 



their composition. During collisions where the relative ve- 
locities exceed 1 m s~^, inorganic compound solids larger 
than 10 cm are destroyed, rather than sticking together and 
becoming a larger, aggregate body (Marov et al. 2008). In 
addition, in circumstellar discs, the time taken for a metre- 
sized solid to fall to the central body is of the order of 100 
years, due to gas drag. Thus, the growth of such solids in the 
disc may take only some tens of rotations (Wiedenschilling 
2000, Armitage 2006). If the growth goes slower, the pro- 
toplanetary disc will lose a large proportion of its 'solid' 
matter, which is the material needed for the formation of 
Earth- group planets, and asteroids and comets. 

The formation of planetesimals can occur either in a 
massive circumstellar disc (age 0.1-1 million years, with the 
mass of the disc comparable to that of the central body) or 
in a medium-mass disc (age exceeding 1 million years, with 
the mass of the disc an order of magnitude less than the 
mass of the central body). 

One possible way to rapidly assemble metre-sized boul- 
ders into a body of planetismal size could be via gravita- 
tional instabilities triggered by the collective motion of the 
solids subdisc. 

Gravitational instability in discs with a central body 
has been investigated since the 1950s (Edgeworth 1949). It 
was found that to experience local Jeans instability, a disc 
with a central body needs lower dispersions in the velocities 
of solids (Gurevich & Lebedinsky, 1950), as well as lower 
gas temperatures (Safronov 1960) or higher-density matter 
(as compared with a motionless medium). Such conditions 
— which are close to those that trigger other instabilities 
in the disc (Fridman 2008) — may lead to the formation of 
rings, spirals, complex wave structures or individual clumps 
(Snytnikov et al. 2004). The formation of clumps in a ro- 
tating medium is hindered by gas heating, as well as the 
increasing velocity dispersion of solids produced by gravita- 
tional field perturbations in spiral and other wave structures. 
The non-linear behaviour of instabilities and the appearance 
of clumps (whose densities increase under self-gravitation) 
are therefore of interest. In any case, for discs to experi- 
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ence fragmentation, they should be quite dense and cold. 
Scenarios for the development of instabilities (Levin 1965, 
Goldreich & Ward 1973) suggest that stable discs transform 
into unstable ones due to the sedimentation of solids on the 
equatorial plane. For instability development it was neces- 
sary to provide the surface density of the primary solids crit- 
ical for this particle velocity dispersion. However, the forma- 
tion of a dense primary solids subdisc is prevented by solids 
scattering, which is caused by turbulent gas flows result- 
ing from hydrodynamic instability of the Kelvin-Helmholtz 
type (Cuzzi, Dobrovolskis & Champney 1993). This instabil- 
ity arises because both components of the two-phase subdisc 
located in the equatorial plane rotate at the Keplerian ve- 
locity, whereas the angular velocity of gas outside the plane 
is lower, due to the radial pressure gradient. A substantial 
difference in subdisc thickness (with respect to the gas and 
condensed phases) is not therefore expected. 

On the other hand, a disc can reach its instability 
threshold not only by increasing its density, but also via 
a decrease in temperature. The main problem with imple- 
menting this mechanism in a medium-mass disc is that the 
lowering of the gas temperature in the disc is accompanied 
by a rapid formation of spiral structures, which do not trans- 
form into clumps. According to estimates made by Makalkin 
& Zighna (2004) and Marov et al. (2008), at the stage of 
planetesimal formation in the Solar System, the temperature 
of the gas component in the disc at the Earth's radius was 
500-700 K when the disc had uniformly distributed solids, 
and 250-500 K after the formation of a thin subdisc of pri- 
mary solids. Spiral waves appear at higher temperatures in 
the medium-mass disc. In these spirals, the dispersion in the 
velocities of solids increases (Rice et al. 2004). These struc- 
tures — which result from the development of gravitational 
instability in the gas — scatter the subdisc of primary solids 
with sizes exceeding 1 m (Rice et al. 2006). Rapid cooling, 
which typically occurs in a time comparable to the time 
taken for the disc rotation, was considered to be a necessary 
condition for clump formation in such a disc (Gammie 2001, 
Rice et al. 2003). This imposes stringent restrictions on the 
mass and temperature of the disc, as well as on the zones 
of formation of such clumps (Rafikov 2009). Further studies 
revealed that such restrictions could not provide conclusive 
evidence, proving — or disproving — that gravitational in- 
stability is the mechanism responsible for the formation of 
large bodies in the discs. The reason is that the critical time 
of cooling depends on the 7 - constant ratio of the spe- 
cific heats (Rice et al. 2005), the history and rate of cooling 
(Clarke, Harper-Clark & Lodato 2007), and the ratio of the 
mass of the disk to the mass of the central body (Meru & 
Bate 2010). Thus, for a disc whose mass is 10 times smaller 
than that of the central body, a cooling time equal to 15 
disc rotations is now taken to be critical (Lodato & Clarke 
2011). However, in the case of a medium-mass disc, such gas 
cooling conditions can hardly be expected to occur, due to 
the growing radiation from the star. 

The mechanisms that facilitate the formation of plan- 
etesimals from large bodies have previously been studied for 
medium-mass discs. According to the computer simulation 
made by Rice et al. (2004, 2006) for a two-phase disc, metre- 
sized solid particles can concentrate in spiral arms for some 
time, due to gas dragging. Computational experiments by 
Youdin & Shu (2002) demonstrated that a turbulent gas 



flow resulting from the difference in the angular velocities of 
the primary solids subdisc and the gas disc did not prevent 
gravitational instability which increases the concentration 
of solids. A condition for the development of instability is 
determined by the compression of the solid phase subdisc to- 
wards the star, in the equatorial plane. Such compression in- 
creases the ratio of the mass densities of the solid component 
and the gas by a factor of 2-10. The density of the medium 
in these regions starts to grow under the action of a self- 
gravitational field. Marov et al. (2008) noted that the pos- 
sibility of a local concentration of bodies should be consid- 
ered; this local concentration may result from the differential 
rotation of the gas subdisc in large turbulent structures. Ac- 
cording to Cuzzi et al.(1993), such bodies will be represented 
by metre-sized boulders, since solids exceeding 1 m in size 
settle most efficiently on the equatorial plane. Meanwhile, 
smaller agglomerates are removed from this plane by a tur- 
bulent fiow that results from the different angular velocities 
of the gas and solid-phase subdiscs. 

Overall, the problem of planetesimal formation from 
primary solids, boulders and agglomerates in medium-mass 
quasi-stationary discs has not yet been convincingly solved. 
We therefore examined the possibility of planetesimal forma- 
tion in the intermediate period between the existence of the 
massive accretion disc (age of more than 0.1 million years, 
with the mass of the disc comparable to the mass of the 
central body) and the medium-mass disc (age of more than 
1 million years, with the mass of the disc approximately 10 
times smaller than the mass of the central body). 

For massive disc formation Machida, Inutsuka & Mat- 
sumoto (2008) simulated the structure of the gas fiow and 
demonstrated that the streams outside the disc plane along 
the rotation axis must exist. These outgoing fiows can pro- 
vide intense non-radiative cooling of gaseous disc under the 
condition of weak radiative heating from protostar of a small 
mass. 

Johansen et al. (2007) investigated local gravitational 
collapses in medium of sub-metre-sized solids in the disc. 
In their simulations the dispersion in the velocities of such 
solids was determined by turbulent fiuctuations in the gas 
fiow. Owing to the primarily collective motion of the gas 
and solids, the dispersion may be as high as ten metres per 
second. Johansen et al. (2007) found that in a local region 
with a self-gravitating medium and an external gravitational 
field exerted by a protostar, the drag force resulted in gravi- 
tational collapse in a subsystem of solids, over several orbital 
periods. Gas was not involved in such collapses. This mech- 
anism was able to explain the growth of sub-metre-sized 
agglomerates into bodies with a size of 10 metres. These 
10-metre-sized bodies experience rare but high-velocity col- 
lisions with each other during orbital time. In contrast to 
sub-metre-sized solids, drag forces do not cause such en- 
larged bodies to move in concert with gases. The collective 
motion of these bodies in the disc is described by the Vlasov 
equation. 

The main thesis of this paper is that the gravitational 
interaction of the two phases (gases and primary solids over 
a metre in size) affects the stability of the entire disc, and 
changes the conditions necessary for the formation of clumps 
in such a system. We investigated the gravitational interac- 
tions between the gas and solid phases in a massive disc, 
during the development of an instability that could lead to 
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the formation of clumps. As it was explained above the sur- 
face density of the solid phase was assumed to be greater 
than 1 per cent of the gas surface density. The solids sub- 
disc reaches its fragmentation threshold only due to presence 
of massive gaseous envelope. 

The second section of the paper describes a mathe- 
matical model for the two-phase protoplanetary disc during 
clump formation, and the third section deals with the nu- 
merical algorithms and code we used in the computational 
experiments. The fourth section presents estimates for the 
effective Jeans lengths of the gravitational interactions be- 
tween the gases and the primary solids subdisc. The fifth 
section shows the results of numerical modelling for the dy- 
namics of the two-phase system with a central body, for time 
periods covering several disc rotations. 



2 A MATHEMATICAL MODEL OF THE 
PROTOPLANETARY DISC DURING 
CLUMP FORMATION 

2.1 Basic equations 

The computational experiments reported in this paper were 
carried out within a quasi-3D model of the disc. This means 
that we neglected the vertical motion of matter, and con- 
sidered the dynamics of the razor-thin disc where its entire 
mass was concentrated inside the equatorial plane of the 
system. The gravitational potential of the multiphase disc 
was calculated using the 3D Laplace equation. 

The limits of such model applied to the disc description 
were discussed in detail by Fridman & Khoruzhii (2003). 
They found that the model is valid, first, when the typical 
length of density perturbation in the disc is greater than its 
thickness; secondly, when the changes of the dynamical pro- 
cesses time are slow with respect to orbital time. The first 
condition is satisfied for all process in the disc except for its 
collapsing. The second condition is satisfied for gravitational 
instability almost everywhere inside the disc except for its 
outer edge where dynamical and orbital time are compara- 
ble. 

The gas component was described by the following gas 
dynamics equations: 



da 

dt 



+ div{av) — 0, 



dt 



+ (,;,V)5'* = 0, p* =TV. 



These gas dynamics equations include surface quantities 
that were obtained from volume quantities by integration 
with respect to the vertical coordinate z: 



C = CTnas = 



/+00 r-\-oo 
Pgasdz; p* = pdz. 
-oo J —oo 



Here, v — (vx^Vy) is the two-component gas velocity, 

* T* 
and p* is the surface gas pressure. T* = — , 5** = In — 

are the quantities similar to gas temperature and entropy. 7* 
is a 2D version of 7(Polyachenko & Fridman 1976), which is 



related to the constant ratio of specific heats as 7* =3 . 

7 

^ - is the gravitational potential in which the motion occurs. 

In our model, the solid-phase subdisc was represented 
by 1-10 m solids moving in such a manner that two solids col- 
lided with a frequency not exceeding one event per rotation 
around the protostar. The dynamics of the primary solids 
subdisc were described by the Vlasov equation, neglecting 
collisions between solids at times longer than several rota- 
tions: 

dt dr du ' 

where a = — V$, a is the acceleration of particles in an 
external and self-consistent field, u = (ur^Ucf)) is the two- 
component velocity of the particles, and / = /(t, r, u) is a 
function of the particle velocity distribution at a point in 
the disc with the coordinate r. This function is related to 
the surface density of the particles as a par — J fdudz. Note 
that the radius of the solids does not appear explicitly in 
the Vlasov equation; however, the description of the primary 
solids subdisc as a collisionless system implies a certain size 
range for these solids. 

$ is the gravitational potential, which is the sum of the 
potential of the motionless central body and the potential 

Mc 

of the razor-thin disc, $ = $1 + $2,^1 — , where Mc 

r 

is the mass of the central body. ^2 is the potential of the 
self-consistent gravitational field, which is defined in general 
from the Dirichlet problem for the 3D Poisson equation 



A$2 = 4:7r{ppar + Pgas 



$2 



\/r 



0. 



For ppar + Pgas = {crpar{r, 0) + crgas{r, (l)))6{z), whcrc S{z) is 
the Dirac delta function, $2 is found as a solution of mixed 
problem for the Laplace equation 



A$2 

d^2 , 
dz ' 



0, <l>2 



0, 



The equations are written using dimensionless variables. 
The basic dimensional quantities are G (the gravitational 
constant), and Ro = 10 au. Mo = Mq = 2 • 10^° kg, which 
are the typical size and mass of the system. The gas compo- 
nent and the boulders subdisc interact through a common 
gravitational field. 

2.2 Initial conditions 

The surface temperature and density of the disc were speci- 
fied at zero time. In the calculations presented in this paper, 
the density of the gas and primary solids subdisc was taken 
as a Mackloren disc of mass Mpar,gas and radius R: 



3Mp, 



27ri?2 



0, R. 



The gas temperature at zero time was specified as 



T (r) ^ {1 — —) OT T (r) ^ a{r), using a given To, which 
R 

is the temperature in the 'centre' of the disc. 

The initial velocities of the solids were specified as the 
sum of regular and random components, u = u'-\-u'\ where 
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is the regular velocity, and u'^ is the random velocity. The 
gas velocity and the regular velocity of particles were de- 
termined from an equilibrium between centrifugal and cen- 

tripetal gravitational forces: = + ^, = 

r a or or r or 

where Vr = 0, and = 0. The random velocity of the 
particles i^" was specified by the Gaussian law with a zero 
average and a prescribed dispersion Vd- 



3 NUMERICAL METHODS AND SETUPS 

We carried out our computer simulation of protoplanetary 
disc dynamics using a Sombrero code based on the method 
of splitting with respect to the physical processes involved 
(Stoyanovskaya & Snytnikov 2010). The Vlasov equation, 
the gas dynamics equations, and the mixed problem for the 
Laplace equation were solved at each time step. 

The solution of the Vlasov equation was obtained using 
the particle-in-cell (PIC) method (Hockney 1987). We com- 
bined it with the grid method solution of the Laplace equa- 
tion. We interpolated the surface density of solid bodies on 
the regular polar grid using the PIC kernel. The equations 
of each particle motions are characteristics of the Vlasov 

dit dl^ 
equation: — — = — V$, — ;— = it. They were integrated 

dt dt 
by leap-frog scheme (Snytnikov et al. 2004). 

The gas dynamics equations were solved using the SPH 

method (Monaghan 1992). The SPH calculation formulas 

implemented in Sombrero were obtained from quasi-3D gas 

dynamics equations, written in the Lagrangian form: 



da 
Itt 



• divv, 



d d 



dv 1^ dr dS 

—- = Vp-V<I>, — =v, — - 

dt a dt dt 



0, 



where — = — v \/. 
dt dt 

We used the cubic spline for 2D space as a kernel W: 



W{q,h) = 



[(2 -^)3_ 4(1-^)3] 
[2-#, 

0, 



^ ^ ^ 1, 

1 ^ g ^ 2, 

g>2, 



where q = — — , x - radius vector of a space point. 
h 

The surface density of the gas where the zth parti- 
cle resides was calculated as the sum interpolant ai = 
^^^iTTijWij, where N was the number of simulated SPH 
particles. We used the adaptive smoothing length hi = 

2 4 / . ^ , which was implemented through arithmetic 



10-3 +cr. 
kernel averaging 



r, -7;(W(\rr-T,\,h^) + W(\n-T,\,h,)). 



The equation of motion was approximated such that the 
impulse and angular momentum were preserved: 



Hz, = 



0, 



Pij 



1 , 



For calculation formulas we used notations: 



Wij = W{\ri 



h 



dq 



We used a standard artificial viscosity with parameters 
a = 1, P = 1 (Monaghan 1992) in Sombrero code to make 
the calculation of shearing motion possible. In our model 
we used the entropy equation, which doesn't describe the 
shock waves. That is why we introduced an artificial viscos- 
ity to the motion equation only. The kinetic energy of gas 
decreased due to the viscosity on 10 per cent 100 x 128 x 100 
grid cells and 40000 SPH particles for the test calculation. 
The contribution of artificial viscosity is decreased by in- 
creasing number of SPH particles. Numerical experiments 
show that the velocity field in clumps is not affected by 
the number of SPH particles. If we change the model by di- 
rect artificial heating it provides significant increasing of gas 
temperature and the transition of its fiow from supersonic 
to subsonic in the inner part of the disc. 

The Laplace equation was solved using a cylindrical re- 
gion of space, with its lower face determining the disc plane. 
The particles simulating the gas disc and the primary solids 
subdisc were located in this plane. The Laplace equation was 
solved for the entire volume of the cylinder, using a cylindri- 
cal coordinate system. A boundary condition fixing the po- 
tential to zero at infinity was transferred to the side surfaces 
of the cylinder and the upper edge of the calculated region at 
each time step, via the decomposition of the potential into 
multifields (up to quadrupole moments). The radius of the 
calculated region was typically twice as large as the initial ra- 
dius of the disc, and the height of the cylinder was chosen to 
be equal to the radius of the calculated region. The Laplace 
equation was solved using an iterative combined method, 
where a value from the previous time step was taken as an 
initial approximation. This method employed fast Fourier 
transforms in angular coordinates, combined with a succes- 
sive block over-relaxation procedure (Snytnikov et al. 2004). 

The gravitational force affecting the SPH and PIC par- 
ticles was calculated by a linear interpolation of the force 
rate in the mesh points which was obtained from the mid- 
plain gravitational potential value. 

Most of the calculations presented below were per- 
formed on a 100 X 128 X 100 grid. The particles did not reach 
the boundary of the calculated region in any of our calcula- 
tions. To verify that the changes observed in the solutions 
were due to the parameters for the low-mass component of 
the disc (and were not related to numerical effects) , all calcu- 
lations were performed using the same numerical algorithm 
parameters. The gas disc was represented by 4 x 10^ SPH 
particles. In this case, the mass of the particles within the 
smoothing length did not locally exceed the corresponding 
Jeans mass (Bate & Burkert 1997). The primary solids sub- 
disc was represented by 5 x 10^ PIC particles. In the disc 
plane, the grid size was [r, 0] = 100 x 128 with 50 x 128 
cells representing the disc at zero time. The average number 
of PIC particles in each cell in the disc plane was there- 
fore 1000, which provided a ^ 3 per cent fiuctuation in the 
density and other calculated quantities. 

The resolution requirements for multi-phase and single- 
component medium studies are similar. Usually the Jeans 
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length in gas and subdisc of solids is greater than a few grid 
cells and smoothing lengthes of SPH particles. In the next 
section we discuss existence of effective Jeans length for gas- 
solid bodies medium. Hereby we suppose that a maximum 
grid cell size and a smoothing length must be less than the 
effective Jeans length. 

In addition, the applicability of the methods imple- 
mented in Sombrero for the simulation of the dynamics of 
axisymmetric and radial-azimuthal structures (that form in 
the two-phase medium of the gravitating disc) was verified; 
this was achieved by comparing the results of computational 
experiments performed using the SPH and FLIC methods 
(Stoyanovskaya & Snytnikov 2010). It was shown that the 
SPH method was able to reproduce nonlinear waves in the 
gas and collisionless solids medium, with shear and counter 
flows emerging in the system. 



4 GRAVITATIONAL INSTABILITY IN A 
MEDIUM OF GAS AND COLLISIONLESS 
PRIMARY SOLIDS MEDIUM: ANALYTICAL 
EXPECTATIONS 

The scale separating the growing and decaying perturba- 
tions in a gravitating medium is determined by the local 
Jeans length. The Jeans length is obtained from an analysis 
of the dispersion ratios, which gives an estimate of the mini- 
mum density perturbation in a system developing under the 
action of a self-gravitational field. The long-wave stability of 
a collisionless, fiat layer was considered by Fridman & Poly- 
achenko (1984). They demonstrated that the critical length 
of the wave is equal to 



The vibrations of the incompressible gas layer (in the same 
long- wave approximation,) are expressed by the relationship 



A - — 



sfJi GCFg 



Similar expressions for the critical extent of the distortions 
are available for the Toomre instability in a heterogeneous 
rotating disc. The dynamics of gases and collisionless pri- 
mary solids under a common gravitational field (in the same 
long-wave approximation of a fiat layer in a circumstellar 
disc, for perturbation of a potential $ near the equatorial 
plane, in scale A) can be written as 4> ~ G{apar + crJas)A, 
where apar, and a gas are perturbations in the particle and 
gas surface densities. For rapidly growing waves (rapid com- 
pared with the rotation period around a protostar), with 
reference to 

1 agas (Jgas^ 2 ^ 

— ~ , VdCFpar ~ CFpar^, 

. . 1 1 1 

one can obtam — ~ h , or 



Agfas 



A A.par Agras 
Ar)ar 



1 + ^^ 



where rUparVd is a typical 'temperature' for the velocity dis- 
tribution function of solids (neglecting motion perpendicular 



to the disc plane), rupar is the mass of a solid particle, and T 
and /i are the temperature and molecular weight of the gas. 
A gives the critical wavelength for a two-component disc, 
the so-called effective Jeans length for a two-phase system. 



Relation Q = 



Aj 



< 1 gives the classical Toomre 



value for the gravitational instability in the rotating disc, 
where A^ot = —pSr^^ ^ - epicyclical frequency. For mul- 
tiphase disc a = CTpar + cr^as, which defines Arot = 
. Then, using the obtained value of effec- 



G{o'gas ~\~ CTpar) 



tive Jeans length the Toomre parameter can be written as 

Q 



par ^gas 

^gas+^par G{(7gas + CTpar) ' 

Inside the disc, for a mixture of hydrogen and he- 
lium at a temperature T ^ 300 K, Cs ~ 1000 m s~^. For 
CTpar /cr gas ~ 5 * 10~^ (a massivc gas disc), and with a dis- 
persion in the solid velocities of Vd ~ 100 m s~^, the length 
A decreased 6- fold compared with Agas] for Vd ^ 20 m , a 
100-fold decrease was observed. It was therefore shown that 
the presence of a second low-mass component could strongly 
decrease the effective Jeans length for the two-phase system 
A. By developing the aforementioned instabilities in its sub- 
system, the low-mass phase was able to facilitate the de- 
velopment of gravitational instability in the entire medium. 
The role of the gas was reduced to providing a drag force 
on individual solids; such friction decreased the velocity dis- 
persion of the primary solids. The gas therefore served as a 
massive medium, in which density perturbations transferred 
from the primary solids subdisc could develop, and produce 
clumps. The physics of the process is similar to the cooling 
of heavy ions by light electrons (which move in concert with 
the heavy ions) that occurs in accelerator physics (Nikitin, 
Snytnikov & Vshivkov 2004). 



5 COMPUTATIONAL EXPERIMENTS 

We will now clarify how the gravitational interactions be- 
tween the cloud of low-mass collisionless primary solids and 
the massive gas can affect the formation and development of 
structures in the self-gravitating disc, on its rotation around 
the central body. Three cases can be distinguished: first, 
structures can form and develop due to instabilities in the 
gas, and can subsequently involve primary solids in their dy- 
namics. Second, structures can be initiated in the primary 
solids subdisc, and can remain there for their entire lifetime, 
without any noticeable effect on the gas. Third, structures 
can be generated in the gas, but under the development of 
'two-phase' instability, where the overall dynamics of the gas 
are affected by the low-mass primary solids subdisc. 

The last mode is the most interesting; we will demon- 
strate that the 'two-phase' system (gas containing a subdisc 
that has low solid velocity dispersions) can act as a source of 
the instability necessary for the formation of self- gravitating 
gas-solid clumps in the disc. 

In order to demonstrate that all three modes exist in 
massive protoplanetary discs, we reproduced the dynamics 
of some rotations of a disc with mass M = 0.55Mo and ra- 
dius R = 2Ro — 20 au, rotating around a central body with 
mass Mc = 0.45Mo. It should be noted that when going 
from a 3D disc to its quasi-3D approximation, we retained a 
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Table 1. Parameters and corresponding Jeans lengths for computational experiments. 



Experiment No 


1 


2 


3 


4 


Mass of gas 


0.52Mo 


0.495Mo 


0.52Mo 


0.52Mo 


Mass of solids 


0.03Mo 


0.055Mo 


0.03Mo 


0.03Mo 


Initial gas temperature at R = 10 au (K) 


610 


1335 


610 


610 


Initial velocity dispersion of solids (m s"-*^) 


1900 


285 


475 


95 


Initial Jeans length in gas 


0.42 


0.97 


0.42 


0.42 


Initial Jeans length in solids at = 1 au 


6.45 


0.2 


0.4 


0.016 


Outcome of instability development 


3 spiral 


10 spiral arms 


5 spiral 


Gas-solid clumps 




arms 


in solids 


arms 


in the 




in gas 


transformed into 


in gas 


inflection points 




and 


3 spiral arms 


and 


of 5 spiral 




solids 


in gas and solids. 


solids 


arms 



Thermalization of solids 



fixed ratio between the masses of the central body and the 
disc, Mc/M. All calculations reported in this paper used 
7* = 5/3. The ratio of the surface densities of the primary 
solids subdisc and the gas disc {(Jpar/crgas) was varied in 
the range of 0.01 - 0.1. This differed from the ratio between 
the solid and gas phases in molecular clouds (which has a 
value between 0.01 and 0.02), due to the accumulation of 
the solid phase in the disc, and the greater thickness of the 
gas subdisc compared with the solids subdisc. In the quasi- 
3D model, (Jpar/o-gas = mpar/rrigas- By specifying the ratio 
of the surface densities and the mass of the entire disc, we 
determined the masses of the gas and the primary solids 
subdisc. Thus, in calculations 1, 3, 4 (Table 1), the mass of 
gas in the disc was Mgas = 0.52Mo, and the mass of the par- 
ticles was Mpar = 0.03Mo, which corresponds to the surface 
density ratio apar/(Jgas ~ 0.058. In contrast, in calculation 
2, Mgas 0.495Mo and Mpar 0.055Mo, which gave a 
surface density ratio of g par I (J gas ~ 0.11. 

We specified the density and temperature (dispersion) 
of the disc to give an initial state of unstable equilibrium. 
Such a choice of initial state made it possible to follow the 
main stages of the development of instability in the calcu- 
lations. The initial gas temperature and solid velocity dis- 
persion conditions were set to provide a certain level of in- 
stability in each component of the disc. Table 1 provides a 
summary of the parameters and corresponding Jeans lengths 
for each run. Fig. 2 (which shows the initial radial distribu- 
tions of the temperature and angular velocity in the gas, for 
experiment 1) shows that ultrasonic gas flow with differen- 
tial rotation was reproduced. 

Experiment 1 was run in a mode where the formation 
of the disc structure was determined by the gas, while the 
gravitational field of the massive component meant that the 
solids were also involved in the structures. In this run, the 
initial Jeans length A for the gas (which was constant over 
the entire disc) was equal to 0.42, which corresponded to a 
gas temperature profile T* ^ Ggas- For such a dependence, 
at a radius R— 10 au the temperature was 610 K; at a radius 
R— 1 au, the temperature was 700 K. The Jeans length in 
the primary solids subdisc increased from the centre to the 
periphery, and reached 6.45 at = 1 au, due to an initial 
primary solids velocity dispersion of Vd — 1900 m s~^. 



800 



600 



400 



200 




Values for the Toomre parameter Q 



CsK 
1^ 



(which 



characterises the level of gravitational stability of a rotating 
flat disc, and where k'^Vl for discs with near-Keplerian ro- 



15 ^ 20 

r(au) 



Figure 2. Gas temperature and angular velocity versus radius at 
initial time, for experiments 1, 3, and 4 



tation) are shown in Fig. 3 for calculations 1, 3 and 4, for 
gases and primary solids. At a specified gas mass and tem- 
perature, a value of Q < 1 corresponded to the disc region 
R > 2 au, which indicated that the radial-ring instability 
could develop in the gaseous disc. 

As shown in Fig. 4 in the images from Run 1 (the im- 
ages show the surface densities of the gas and particles), 
broad density rings had formed in the gaseous disc and on 
its periphery by the time t = 10. The density distribution of 
solids was similar to that in the gas, except for the narrow 
solid density rings at a radius of ~ 8 au. By the time t = 15, 
azimuthal instability had developed, resulting in the forma- 
tion of a three-arm spiral structure in the gas. The width 
of the solid spirals was much smaller than the width of the 
gas spirals. The solids concentrated in the spiral gas arms 
formed their own structure, which was thinner than the gas 
structure. At a later time of t = 23, this relatively thin solid 
structure disappeared, leaving a well developed three- arm 
structure; in this case, the dynamics of the disc were de- 
termined by the gas. Therefore, for high initial velocity dis- 
persions, the bodies showed neither large inhomogeneities in 
their density nor any significant thermalisation of velocities, 
leading to a smooth spiral structure that covered the solids; 
this was the result when the disc dynamics were determined 
almost completely by the massive gas component. 
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Figure 4. : The formation of spiral arms due to the development of gravitational instability in the gas. The surface density of the gas 
(top) and the primary solids subdisc in experiment 1, for time points T = 10, 15, and 23. 

spiral arms. Radial-ring structures in the solid disc were re- 
vealed in the gas. By t = 20, further interactions between 
the spiral waves in the solid phase had produced a distinct 
ten-arm structure in the outer part of the subdisc, and a 
complex wave structure in the interior of the subdisc. Such 
structures were generally absent in the gas component. The 
ON 

function AN = -- — Avr, where A/'-the velocity distribution 

OVr 

function of solids, (see Fig. 6) remained almost unchanged 
for t = 20. By t = 30 (which roughly corresponded to the 
next disc rotation over the periphery), the solids were scat- 
tered in patterns determined by their density waves, with 
their distribution function showing an increase of nearly 
an order of magnitude. Such a distribution function cor- 
responded to a 3- to 4-arm density of solids in the disc (see 
Fig. 5) and weakly revealed structures in the gas. The num- 
ber of particles with significantly lower velocities (obeying 
the distribution function dv < —0.4) was still increasing, 
as was the number of solids with slightly higher velocities 
(with c/t' < 0.2); by t = 50, the disc composed of both solids 
and gas attained an equilibrium thermalised state. Thus, for 
this run, the thermalisation time of the disc (whose insta- 
bility was higher in the solid component than in the gas) 
was smaller than the time taken for the solid subdisc to de- 
cay into separate clumps, under Jeans instability working 
against a moving but uninvolved gas. 

Experiments 3 and 4 showed the development of two- 
phase gravitational instability, with both the massive gas 
component of the disc and the primary solids subdisc in- 
fluencing the formation of structures. The disc dynamics 
calculations were performed using gas parameters similar to 




0.0 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 \- 

Figure 3. Initial value of the Toomre Q parameter versus ra- 
dius, for experiments 1, 3, and 4. 1, 2, and 3 show curves for 
the primary solids subdisc: 1 - experiment 1 (corresponding to 
an initial thermal velocity = 1900 m s"-*^), 2 - experiment 3 
{vd = 475 m s"-*^), 3 - experiment 4 (vd = 95 m s"-*^). 4 - gas 
curve corresponding to all three experiments. 



The benchmark Run 2 — where the gas temperature 
increased more than twofold, and the solid velocity disper- 
sion decreased by a factor of almost 7 compared with Run 1 
— showed radically different disc dynamics. This run corre- 
sponded to the second case considered above for the interac- 
tion between the particles and gas. For such disc parameters, 
the effective Jeans length A was determined by the primary 
solids subdisc, as illustrated in Fig. 5. By t = 10, a strong 
azimuthal instability had developed in the low-mass com- 
ponent of the disc, resulting in the formation of numerous 
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T=10 




T=20 



r j f^ T=30 



T=50 




Figure 5. The formation of structures and the development of gravitational instability in a mode where the value of the Toomre 
parameter for the primary solids subdisc was lower than that for the gas. Logarithm of the primary solids subdisc surface density, for 
time points T = 10, 20, 30, and 50, for experiment 2. 
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4x10 



3x10 
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Figure 6. Function AA^ : 



dN 

dVr 



Avr, where N-the velocity distri- 



bution function of PIC particles in a solid subdisc residing in a 
ring [rl, r2] = [0.45, 0.9], for time points T = Avr = 0.0033, 
T = 20 Avr = 0.007, T = 30 Avr = 0.011, T = 50 Avr = 0.013, 
in experiment 2.) 



those used in experiment 1; the only change was that the 
primary solids velocity dispersion was decreased to 475 m 
s~^ in experiment 3, and 95 m s~^ in experiment 4. 

In Fig. 7, for Run 3, the density distribution in the gas- 



solid subdisc indicated that the development of this process 
in a gas with a colder low-mass component produced a 5-arm 
structure instead of a 3-arm one, after the same amount of 
time. The lateral linear dimensions of the individual density 
waves in the solids spiral (i.e., the thickness of these enti- 
ties) were much smaller than those obtained in experiment 
1 for the spiral waves formed mainly by the gas component. 
Such a decrease in the thickness of the solid 'filaments' also 
took place in the radial-ring structure on the disc periphery, 
which was less clear for the gas. The peripheral structure 
of the solids was characterised by a strong interaction with 
the spiral waves in the middle of the disc, which led to the 
breakdown of the azimuthal symmetry. This interaction pro- 
duced high-density zones, due to the coupling between the 
particles and gas in the radial-ring and spiral- wave domains. 
This result raises the possibility that self- gravitation suffi- 
cient to oppose thermalisation factors in the two-phase disc 
may be triggered in such zones. 

In order to assess this possibility, let us consider the 
results of experiment 4, where the solid velocity dispersion 
was decreased again by a factor of 5, to approximately 0.1 
Cs. In this case, the spiral structure did not change, and the 
spiral retained its full number of arms. This confirmed the 
importance of the role of the massive and unstable gas disc 
in the appearance of the spiral structure. However, over a 
narrow range, the number of arms depended on the solid 
subdisc parameters. The thickness of the solid-phase arms 
decreased significantly in comparison with experiment 3. In 
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Figure 7. The surface density of gas (top) and the logarithm of the surface density of the primary sohds subdisc, for experiments 1, 3, 
and 4, for time points T = 20, 19, and 22, respectively. 



the radial-ring and spiral- wave coupling region, more clumps 
were formed with high densities of both solid and gas. Here, 
the solid phase density exceeded the background values by 
a factor of at least ten. The clumps formed in the regions of 
minimum Jeans length. Thus, the development of instabili- 
ties in the two-phase disc resulted in the formation of five ar- 
eas of increased gas-solid density, with a possible subsequent 
local gravitational collapse of gas and solids. There was in- 
sufficient time for the scattering of solids by density waves 
to occur in these regions; the time was also not sufficient 
for the primary solids subdisc to experience complete ther- 
malisation, or to pass out of its quasi-equilibrium state. This 
conclusion neglects the drag force between gas and solids, an 
assumption that is valid for some tens of rotations. At long 
times, drag due to the gas should be taken into account, as 
it can decrease the velocity dispersion of solids at a certain 
radius from the central body. 

To make sure that fragmentation in experiment 4 is 
not a numerical artefact we calculated this regime with in- 
creasing numerical resolution. We consequently increased 
the number of SPH and PIC particles four times and the 
number of meshpoints twice in each direction. Fig. 8 shows 
the surface density of gas and solid bodies in linear and log- 
arithmic scale respectively. The spiral structure is formed 
almost simultaneously in calculations with increasing nu- 
merical resolution. We calculated the development of phys- 
ical instability. Generally it shows a discontinuous depen- 
dency of solution on problem parameters (Vshivkov, Nikitin 
& Snytnikov, 2003). So details in density plots are changed 
with the numerical resolution variation, when the picture al- 
most remains. All plots demonstrate the appearance of gas 



clumps bound with solid clumps. Increasing the number of 
PIC and SPH particles and meshpoints allowed us to repro- 
duce the wave interaction of the structures with the groups 
of particles. 

These results can be used to explain the mechanism by 
which gas giants are formed. The formation of gas clusters 
due to gravitational instability has been considered to be 
a possible mechanism for the formation of Jupiter and Sat- 
urn (Boss 2000). However, the separate fragmentation of the 
gas disc was feasible (in computational experiments) only 
under strong cooling of the system during several rotations 
(Boss 2000, Durisen et al. 2007, Rice et al. 2006, Meru k 
Bate 2010). Our calculations demonstrated that in a two- 
phase gravitating medium, the low-mass part of the disc 
(composed of metre-sized solids) could initiate and acceler- 
ate the gravitational growth of waves in both components 
of the system, for the case of low solid velocity dispersions. 



6 CONCLUSIONS 

In this study, we examined a major bottleneck in the under- 
standing of the process of planet formation; namely, the for- 
mation of large bodies (planetesimals and planet embryos) 
from metre-sized boulders in the circumstellar disc. A plan- 
etesimal in its emergent stage is considered as a clump of gas 
and solids whose gravitational field preserves its mass when 
the clump moves. We have proposed a mechanism that can 
explain the formation of planetesimals in massive accreting 
discs: in this emergent stage, the mass of a protostar is al- 
most equal to the mass of the circumstellar disc. Gas leaves 
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Figure 8. The surface density of gas (top) and the logarithm of the surface density of the primary sohds subdisc, for experiment 4 
with varying numerical resolution, for time points T = 22, 23, and 23, respectively. The left column - 40000 SPH particles, 2500000 PIC 
particles, 100 x 128 x 100 grid cells, the middle column - 160000 SPH particles, 10000000 PIC particles, 200 x 256 x 200 grid cells, the 
right column - 640000 SPH particles, 40000000 PIC particles, 400 x 512 x 400 grid cells. 



the solid subdisc, being reflected from the equatorial plane, 
and is scattered to outer space. The gas temperature drops 
during this process, under conditions of low radiation from 
the low-mass protostar. Solid bodies grow to metric size with 
the occurrence of collisions, and lose relative velocity due to 
gas drag. In general, the solid component stays in the equa- 
torial plane of the disc, while the gas leaves it. The fraction 
of the solid component therefore increases with respect to 
the gas. The increasing solid density, in combination with 
the decreasing gas temperature, causes the two-phase disc 
to transfer into a marginal state, allowing the development 
of gravitational instability of some type. 

Our computer simulation showed that self- gravitating 
clumps were formed in a massive disc via the development 
of a 'two-phase' Jeans instability in the gas-primary bod- 
ies medium. These bodies were of larger than metre-sized, 
and rotated around the protostar without the occurrence of 
collisions per orbital time. In these unstable conditions, the 
overall gas dynamics were affected by the primary solids 
subdisc, via its gravitational field. This implied that the 
possibility of clump formation was determined both by the 
rate of gas cooling and its density redistribution, and by 
the rate of concentration of large (over 1 m) primary solids, 
and the decrease in their velocity dispersions (cooling of pri- 
mary solids). We found that a velocity dispersion of 0.1 Cs in 
the boulder subdisc was sufficient to cause the formation of 
self- gravitating clumps of gas and boulders. In such regimes, 
the time taken for the formation of collapsing objects was 



less than time taken for boulders to dissipate in the density 
waves of the medium. 
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